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ABSTRACT 

The global development of magnetohydrodynamic turbulence in an accretion disk is studied within 
a simplified disk model that omits vertical stratification. Starting with a weak vertical seed field, a 
saturated state is obtained after a few tens of orbits in which the energy in the predominantly toroidal 
magnetic field is still subthermal. The efficiency of angular momentum transport, parameterized by the 
Shakura-Sunyaev a parameter, is of the order of 10~^. The dominant contribution to a comes from 
magnetic stresses, which are enhanced by the presence of weak net vertical fields. The power spectra of 
the magnetic fields are flat or decline only slowly towards the largest scales accessible in the calculation, 
suggesting that the viscosity arising from MHD turbulence may not be a locally determined quantity. 
I discuss how these results compare with observationally inferred values of a, and possible implications 
for models of jet formation. 

Subject headings: accretion, accretion disks — hydrodynamics — instabilities — magnetic fields — 
MHD ~ turbulence 

For visualization of simulations see http://www.cita.utoronto.ca/~armitage/global_abs.litml 



1. INTRODUCTION 

The Balbus-Hawley instability is the most generally ap- 
plicable mechanism known to initiate turbulence and out- 
ward angular momentum transport in accretion disks (Bal- 
bus & Hawley 1991). This is a linear, local instability that 
exists for rotating flows threaded by a weak magnetic field 
with dO^/dr < 0, conditions satisfied in disks (for earlier 
discussions see Velikhov 1959, Chandrasekhar 1961). A 
vigorous growth rate is obtained for a wide variety of ini- 
tial magnetic field configurations (Balbus & Hawley 1992; 
Ogilvie & Pringle 1996; Terquem & Papaloizou 1996), im- 
plying that the instability is inescapable for ionized disks 
where the field is well-coupled to the gas. 

Extensive numerical simulations have explored the non- 
linear development of the instability within the local, 
shearing box approximation (for a review, see e.g. Gam- 
mie 1998). Such simulations have convincingly established 
that the nonlinear development of the Balbus-Hawley in- 
stability leads to sustained turbulence and significant an- 
gular momentum transport, typically finding a Shakura- 
Sunyaev (1973) a « 10^^ (Hawley, Gammie & Balbus 
1995, 1996; Stone et al 1996; Brandenburg et al. 1995). 
There is some evidence for cyclic behavior that might have 
important implications for disk variability (Brandenburg 
et al. 1996). Equally important has been the final elimi- 
nation of convection (Stone & Balbus 1996), and the near- 
elimination of nonlinear hydrodynamic turbulence (Bal- 
bus, Hawley & Stone 1996), as plausible rival mecha- 
nisms for angular momentum transport in accretion disks. 
Progress has also been made in trying to understand how 
the rich phenomenology of accretion disk variability can 
arise within a dynamo driven disk model (Armitage, Livio 
& Pringle 1996; Gammie & Menou 1998), although much 
more remains to be done in this area. 

There are many further questions that one may hope 
simulations will address, and not all of them are amenable 



to a local treatment. Most obviously, is the angular mo- 
mentum transport in a disk locally determined? What is 
the structure of the spatial and time variability of the disk 
fields, and are they suitable for launching a magnetically 
driven disk wind or jet? Unsurprisingly, the global cal- 
culations needed to investigate these issues are extremely 
demanding, both as a consequence of the larger computa- 
tional domain and, especially, because of the need to sim- 
ulate regions of low density where the high Alfven speed 
severely limits the timestep of explicit numerical codes. 

In this paper, results are presented from a vertically un- 
stratified global simulation of accretion disk turbulence. 
Such a calculation is evidently missing essential physics. 
There is no buoyancy, no possibility of Parker instability 
(Parker 1979), and no magnetically dominated disk corona 
- all features that are expected to arise in a full disk model 
and which may be crucial for the disk dynamo problem 
(Tout & Pringle 1992). However the lesser computational 
demands permit a preliminary investigation of some of the 
important questions raised by previous, local, simulations. 

2. NUMERICAL SIMULATION 

The equations of ideal magnetohydrodynamics (MHD) 
are solved using the ZEUS-3D code developed by the Lab- 
oratory for Computational Astrophysics (Clarke, Norman 
& Fiedler 1994; Stone & Norman 1992a, 1992b). ZEUS 
is a time explicit eulerian finite difference code that uses 
the method of characteristics (MoC) - constrained trans- 
port scheme to evolve the magnetic fields (Hawley & Stone 
1995; Stone & Norman 1992b). For this simulation an 
isothermal equation of state P = pc^ replaces the internal 
energy equation, so the remaining equations are: 
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where the symbols have their conventional meanings. 
There is no explicit resistivity in the calculation, reconnec- 
tion occurs numerically on the grid scale. We use second 
order interpolation (van Leer 1977) for all advccted quan- 
tities, and the latest (and allegedly most stable) version of 
the MoC algorithm for the induction equation. 

2.1. Initial and boundary conditions 

The calculation is performed in cylindrical polar geome- 
try (r, (f>, z), in a volume bounded by Tin = 1, Tout = 4, and 
z = ±H — 0.4. We take $ = ^(r) = 1/r, implying no vor- 
tical component of gravity. The boundary conditions arc 
periodic in z, reflecting at r = rin (fr = Br =0), and set 
to outflow at r = rout- Outflow boundary conditions are 
implemented as a simple extrapolation of fluid variables on 
the grid into the boundary zones. These boundary condi- 
tions admit the development of net vertical flux through 
the computational volume as material leaves the grid - the 
import of this for the resulting a will be discussed later. 
We use a grid of (n^ = 128, = 320, = 32) zones, with 
uniform zoning. 

An initial state is obtained by evolving a disk with con- 
stant surface density and Keplerian rotation in two dimen- 
sions until transients due to pressure gradients and the in- 
fluence of the inner boundary condition die out. This re- 
quires a time = 100, where time is measured in units of 
the orbital period at the inner edge. The resulting equilib- 
rium state has E = constant between i? = 1.4 and i? = 3, 
tapering to a low value at the boundaries. The velocity 
profile is Keplerian to better than 10% over the entire ra- 
dial range. No evidence is found for purely hydrodynamic 
disk instabilities. We then add a weak vertical magnetic 
field that is non-zero only between r = 1.5 and r = 3.5, of 
the form, 



B^ir) = — sin(7r(r- 1.5)). 
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chosen simply to be divergence free and to have zero net 
flux in the z direction. The simulation is then evolved for 
a further 80 orbits. The computational cost of this using 
ZEUS is approximately 2 x 10^^ floating point operations 
- a large but not outrageous number for current worksta- 
tions. 

2.2. Results 

The initial state is immediately unstable, leading to 
rapid growth in the magnetic and perturbed velocity fields. 
Around 40 inner orbits of evolution were required before 
a saturated state was obtained, which was then followed 
for an additional time At = 40 without any further qual- 
itative changes in disk behavior occurring. The magnetic 
energy in the saturated state is dominated by the toroidal 
field component and is subthermal, ^ 0.2 — 0.3 of the ther- 
mal energy, the corresponding energies in Br and B^ are 
respectively 10~^ and 6 x 10~^ of the thermal energy. The 
velocity perturbations are mildly anisotropic, Vz and Sv^j, 
display roughly gaussian distributions with width « 0.3cs, 
while that of Vr has a width of w 0.4 — 0.5cs. The to- 
tal energy in the perturbed velocity field is an order of 
magnitude below that of the magnetic fields. 



Fig. 1 shows the surface density and vertically averaged 
Bz component of the magnetic field at the conclusion of 
the calculation. We plot the overdensity ST, relative to the 
azimuthally averaged surface density profile. 
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The magnetic field is plotted without any such normaliza- 
tion. 

The combination of turbulence and shear leads to a 
ragged spiral pattern in the surface density, though the 
azimuthal fluctuations are relatively small and described 
by a gaussian with a width ~ 0.35E(r). The magnetic 
field displays a filamentary structure, both in projection 
and when visualized in three dimensions. Visually it is 
evident that the field shows considerable coherence over 
large scales, especially in azimuth. Br and B^ (not plot- 
ted) show similar patterns. 

2.3. Efficiency of angular momentum transport 

The efficiency of angular momentum transport in the 
simulation can be measured by a, the ratio of the shear 
stress Wr^ to the gas pressure P in the form. 
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This is consistent with a relation between the shear vis- 
cosity v and the disk sound speed Cg of the usual form, 
u = ac^/fi, with fl the Keplerian angular velocity. The 
magnetic and fluid contributions to the total stress are 
then given by (e.g. Gammie 1998), 
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respectively, where the brackets denote an average in the 
spatial co-ordinates. 

Fig. 2 shows the time evolution of aM and aji at a sin- 
gle radius r = 0.5(rin + rout) in the disk. Apart from 
some transient waves at early times, a from both fluid 
and magnetic stresses is positive, with the bulk of the an- 
gular momentum transport being provided by magnetic 
stresses. At the end of the calculation, the total a at this 
radial location is in the range a = 0.05 — 0.1. 

Evaluating um and an as functions of radius at t = 80, 
an is found to be roughly constant across the grid in the 
range q:_r = 0.01 — 0.02. aM varies by factors of a few, 
and is close to its minimum value at the location in the 
center of the grid used for plotting the time series in Fig. 2. 
Taking an average over the entire simulation volume, we 
obtain, 

aM=i0.17, Qfi~ 2.0x10-2. (8) 

The value of a can also be estimated by comparing the 
simulation to the results of one dimensional diffusive disk 
models (e.g. Pringle 1981). Comparing the evolution of 
the radial center of mass between the two calculations one 
again obtains an a value of ~ 0.1, though as mentioned 
already a is unsurprisingly not a constant in the simula- 
tion. 
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2.4. Influence of net vertical magnetic field 

The choice of outflow boundary conditions at rout was 
motivated by the desire to reduce the ampUtude of re- 
flected waves at the edge of the grid. However these bound- 
ary conditions also aflow the development of net vertical 
magnetic fields that can have a strong influence on the 
properties of disk turbulence. Previous work has found 
that a net field significantly enhances the strength of tur- 
bulence and boosts the derived value of a (Hawley, Gam- 
mie & Balbus 1995). Adopting the parameterization of 
the results of local simulations given by Gammie (1998), 

^,^0.01 + 4^^ + 1^^, (9) 

Cg 4 Cg 

where vaz is the Alfven velocity corresponding to the net 
vertical field, we find that the enhancement of a in this 
simulation due to net fields is expected to be ^ 0.05. This 
is obviously a very crude estimate, as we are using the 
above relation in an untested regime, but it does imply 
that it is premature to conclude that global simulations 
lead to higher values of a than local calculations. 

2.5. Coherence of the m,agnetic field 

The scale of the magnetic fields generated in the disk is 
analyzed via the azimuthal Fourier decomposition, 

CmAr, ^) = ^ fj Bie-'"'Hcj). (10) 

Fig. 3 plots m|CmP for the three field components, in each 
case averaged over the the entire volume of the simulation. 

The azimuthal power spectra show similar shapes for 
each of the magnetic field components. There is a power- 
law decline somewhat steeper than that expected from 
Kolmogorov turbulence at large m, and a break at m ~ 
10 — 20, corresponding to a physical scale of ^ H in the 
center of the grid. This confirms the visual impression 
that the magnetic fields are patchy with typical scales of 
the order of H . However there is also considerable power 
at the largest scales, with the power per logarithmic inter- 
val in m declining only slowly towards low m for all field 
components in the range m = 1 — 10. Similar conclusions 
follow from the radial power spectra. 

3. DISCUSSION 

In this paper, we have reported on a global simulation 
of an unstratified magnetized accretion disk. As expected 
from analytic considerations (Gurry & Pudritz 1994, 1995, 
1996) and local simulations, MHD instabilities generate 
sustained turbulence that leads to outward transport of 
angular momentum. The generated fields possess consid- 
erable power in azimuthal modes of low to, which corre- 
spond to physical scales considerably in excess of -ff , the 
disk semi-thickness, and display a ragged spiral structure. 
The current simulation does not admit the development of 
the Parker instability, which might depress the power on 
large scales, but with this caveat the results suggest that a 
viscosity originating from MHD turbulence may not be a 
locally determined quantity. The dominant magnetic field 
component is toroidal, and the interaction of this fluctu- 
ating internal fleld with a magnetosphere is likely to be an 



important complication to the already complex picture of 
star-disk interaction in magnetic systems (Miller & Stone 
1997; Torkelsson 1998). 

The efficiency of angular momentum transport seen in 
the calculation, parameterized by the Shakura-Sunyaev a 
prescription, is a « 10^^. This is larger than the value 
obtained from local calculations (Gammie 1998) with zero 
net vertical magnetic flelds, though we have noted that 
the influence of a vertical fleld on the current simulation 
is likely to have boosted the value of a significantly. More 
realistic simulations with demonstrated numerical conver- 
gence are evidently required. However there is no obvious 
discrepancy with the values of a inferred for dwarf novae, 
where modeling of disk outbursts suggests a = 0.1 — 0.3 
(Gannizzo 1993), and for Active Galactic Nuclei, where 
the admittedly poorer observations are consistent with an 
a of 10~^ (Siemiginowska & Czerny 1989). Gonversely it is 
hard to see why a viscosity derived from MHD turbulence 
should be two or three orders of magnitude lower in the 
ionized inner regions of protostcllar disks, as required to 
match the timescales of FU Orionis outbursts within disk 
instability models (Bell & Lin 1994). This may call into 
question the self-regulated aspect of the thermal disk in- 
stability picture for FU Orionis events. The generation of 
magnetic fields of large scale may additionally be impor- 
tant for models of jet formation (eg. Blandford & Payne 
1982; Guyed, Pudritz & Stone 1997; Matsumoto & Shi- 
bata 1997; Konigl 1997), which if generated via a disk dy- 
namo would be expected to be most efficient in relatively 
thick disks or advection dominated fiows (Narayan & Yi 
1995). Observations of which systems produce jets appear 
to be broadly consistent with a model in which {H/R) is 
a controlling parameter, though many other possibilities 
are also viable (Livio 1997). 

I thank Jim Stone and James Murray for valuable dis- 
cussions at the start of this work, and the referee for a very 
prompt and helpful report. I am grateful to Mark Bartelt 
for ensuring the availability of the required computing re- 
sources. 
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Fig. 1. — Maps of the disk surface overdensity (left image) and integrated vertical component of the magnetic field, Bz (right image), at 
t = 80 orbits. The surface density is plotted relative to the azimuthally averaged value, ic S(r, 0)/S(r). Typical azimuthal fluctuations in S 
are at the tens of percent level, typical Bz fields are of the order of 10""^ of the thermal energy. The disk rotates clockwise. 
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Fic;. 2. — Derived Shakura-Sunyaev a parameter at the center of the grid. The solid line shows the total a, the dashed and dotted lines the 
contribution from magnetic and fluid stresses respectively. 




Fic;. 3. — Fourier decomposition of the azimuthal structure of the disk magnetic field, averaged over the simulation volume. Prom top 
downwards, the curves show the power spectra for B^, Br and Bz respectively. 



